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I. 


INTRODUCTION 


As commercial interest in space launch vehicles has 
increased over the last twenty years, there has been an 
increased emphasis on improving the efficiency of vehicles 
used for space launch services. This emphasis will continue 
to receive attention as competition forces the providers of 
launch services to decrease costs. 

This thesis is an initial step in the development of an 
exergy methodology tool for the conceptual, preliminary, and 
detailed design of aerospace vehicles. 

In order to determine the efficiency of a particular 
launch system, it is necessary to identify the various 
sources of efficiency loss and the quantity of efficiency 
loss associated with each source within the system. Although 
energy cannot be created or destroyed, it can be transformed 
into unavailable forms. This thesis will present an 
investigation of one method that may be used to examine the 
efficiency loss from two specific areas of a typical launch 
vehicle. This method uses the concept of available energy or 
exergy as the basis for determining the decrease in the 
energy state that occurs in the shock wave and the boundary 
layer of a supersonic vehicle. 

Typical performance analyzes of energy systems use only 
the first law of thermodynamics. This method of analysis is 
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based on the principle of conservation of energy, and 
suggests that an energy balance accounting for all energy 
flow into or out of the system will provide a measure of the 
system efficiency. While the first law approach is extremely 
important in the performance analysis of energy systems, it 
provides only the overall performance or the actual gross 
performance of the system or process. To establish the 
magnitude, location, and type of loss within a system it is 
necessary to use the Second law of Thermodynamics. The 
second law accounts for the irreversibility of all real 
processes and thus provides a method of calculating the 
losses associated with entropy production that occurs. 

The exergy method of analysis uses both the First and 
Second laws of Thermodynamics to determine the energy that 
is available at various points in a system. While the exergy 
method of analysis has been applied to conventional closed 
cycle thermal systems such as power generation and 
refrigeration systems, it has not, to the author's 
knowledge, been applied to the open cycle space launch 
vehicle. 

"Exergy is defined as the work that is available in a 
gas, fluid, or mass as a result of its nonequilibrium 
condition relative to some reference condition."[Ref. l:p. 
33] Since exergy is defined as the available work in a gas, 
fluid, or mass, it is a special case of defining available 
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energy. Exergy is an explicit property, and therefore can be 
calculated at any point in an energy system. 

The exergy method permits the examination of any process 
in relation to the theoretically most efficient manner that 
the process could be carried out within the environment. An 
exergy analysis allows one to quantify the loss of 
efficiency in a process due to the loss of the quality of 
the energy. 

Exergy is calculated relative to a reference condition 
by the following general equation:[Ref. l:p. 34] 

Exergy = <u-u 0 ) - T 0 (s-s 0 ) + lS(v-v a ) - -YL * 

C 

< z ~ z o) yr~T + E + (!) 

Qc'* o 

E i A i F i (3T*-T*-4T 0 T 3 ) + ... 

where the subscript o denotes the reference condition. The 
terms on the right side of the equation are internal energy, 
change in heat energy, flow work, kinetic energy per unit 
mass, potential energy, chemical energy, and radiation 
energy emission respectively. In the analysis of specific 
systems or processes one or more of the above terms may be 
omitted, and only those relevant to the problem being 
considered need be included.[Ref l:pp.34-36] 

This thesis focuses on the quantity of exergy decrease as 
a result of entropy production across the shock wave and as 
a result of boundary layer formation. As can be seen from 
equation 1, an increase in entropy results in a decrease in 
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exergy. The Second law of Thermodynamics tells us that all 
real process are irreversible and that entropy is increased 
in an irreversible process. 
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II. EXERGY DECREASE ACROSS THE SHOCK WAVE 


Supersonic flow over a blunt body will result in the bow 
shock wave being detached from the body as shown in Figure 
1. The shape of the detached shock wave, its detachment 
distance, and the flowfield between the body and the shock 
are determined by the shape of the body and the free stream 
Mach number. [Ref.2:p.128-129] 



Figure 1: Supersonic Blunt Body in Flowfield. 
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Since the velocity of the body is supersonic, the flow ahead 
of the body is unaware of the presence of the body until the 
shock wave is encountered. The shock wave directly in front 
of the body is approximately normal to the free stream and 
consequently the flow behind the shock wave at this point is 
subsonic. The subsonic flow region is defined by the shock 
wave, the body, and the sonic line as shown in Figure l. 

Because the flowfield between a blunt body and its 
detached shock wave contains both subsonic and supersonic 
flow, an exact determination of the flowfield dynamics is 
extremely complicated to formulate. The steady state 
flowfield in the subsonic region is described by elliptic 
equations ard in the supersonic region by hyperbolic 
equations. The sonic line divides these two regions. It is 
the problem of predicting the flowfield in the transition 
region between the subsonic and supersonic flow that makes 
the flowfield calculations for a supersonic blunt body so 
difficult. 

The flowfield behind a detached shock wave cannot be 
precisely determined without the application of 
sophisticated numerical methods. In keeping with the 
development of a design tool, a simplified method for 
determining the approximate shape of the detached shock wave 
that was developed by Moeckel [Ref. 3] is used. This method 
was subsequently examined by Love [Ref. 4] and with some 


6 









modification was found to be valid for blunt bodies 
traveling at supersonic and hypersonic speeds. 

Love introduced two major modifications to Moeckel's 
method. The first of these was in the method used to 
calculate the shock detachment distance from the body. 
Moeckel used a continuity relation to calculate detachment 
distance and shock shape. Love calculates shock detachment 
distance by the following equation: 


= 0.5 C cot $ det 


( 2 ) 


where x'/d' is the nondimensional detachment distance, C is 
a factor that depends on type of body (spherical, 
cylindrical, etc.), and 5 det is the required deflection 
angle for shock wave detachment, from a conical body, for 
the given Mach number (see Figure .2) . 

The second modification made by Love was in the 
calculation of the angle eta (ij) that the sonic line makes 
with the freestream. Love conducted wind tunnel tests and 
compared his experimental results with his theoretical 
calculations [Ref.4:pp.13-15]. For this analysis, the angle 
r] corresponding to l.lsMslO was obtained by developing the 
followirg analytical expression for the values of r\ vs Mach 
number as set forth by Love [Ref.4:p.45]: 

, 3.5901 + 56.9298 *log(M) 

1 ~ 57.29578 
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Figure 2: Representation of Flow With Detached Shock 
and Notation Used in Analysis. 

A. BODY SHAPE 

The shape of the deflecting body is one of the factors 
determining the shock wave shape. For this study, the body 
was patterned somewhat loosely after the Titan III-A launch 
vehicle. This body shape is essentially a blunt cone and 
cylinder arrangement. The dimensions of the launch vehicle 
used in the analysis are shown in Table 1. 


TABLE 1: LAUNCH VEHICLE DIMENSIONS 


Length Overall 

120 Ft. 

Diameter 

10 Ft. 

Nose Radius 

1.5 Ft. 

Nose Cone Length 

10 Ft. 
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B. MACH NUMBER 


Besides the shape of the body, the other factor in 

determining the shape of the shock wave is the freestream 

Mach number. This study determined the shock shape based on 
l.lsMslO. Although this is far short of the Mach number 
required to reach orbit, the purpose of the study was to 

demonstrate the validity of the method, which can be 

accomplished in the listed Mach range. 

C. ASSUMPTIONS 

The method of approximating the shock shape and location 
developed by Moeckel is based on two assumptions: 

(1) The form of the shock between its foremost point 
and its sonic point is adequately represented by 
a hyperbola asymptotic to the free stream Mach 
lines; and 

(2) the sonic line between the shock and body is 
straight and inclined at an angle that depends 
only on the freestream Mach number. [Ref.3:p.i] 

The first assumption states that the shock shape between its 

vertex, or foremost point, and its sonic point may be 

represented by a hyperbola. For this study, the assumption 

was made that this approximation for the shock shape is 

valid past the sonic point of the shock wave. For 

calculation purposes, the shock wave was taken from the 

vertex of the wave to the point where the angle between the 

shock wave and the asymptote that defines the hyperbola 

differ by no more than one per cent. The equation for the 

hyperbola used in these assumptions is 


9 






( 4 ) 



where x is the longitudinal axis of the body, y is the axis 
perpendicular to the x axis, 0 is the cotangent of the Mach 
angle, and x Q is the distance from the vertex of the shock 
wave to the intersection of its asymptotes.[Ref.3:p.1] 

The second assumption made by Moeckel concerning the sonic 
line was modified by Love as previously discussed and was 
used in this analyzes. 

Air is the fluid used for all of the analysis presented 
in this thesis. However, the analysis method is equally 
valid for any gas, and would only require using the 
appropriate value for the ratio of specific heats (y) and 
the new specific gas constant. It has been assumed that the 
fluid behaves as a perfect gas in all of the calculations. 
This assumption is valid for velocities up to approximately 
M=3 to M=5 [Ref 5:pp.17-20]. Above these velocities 
dissociation of the fluid begins to occur due to atmospheric 
heating. If the gas continues to be heated, ionization will 
occur. As the degree of dissociation and ionization 
increases, the perfect gas assumption becomes less valid. 

The calculation of total entropy across the shock wave 
requires a determination of the mass flow rate of air 
through the shock wave. The air mass flow rate calculation 
is based on the density, cross sectional area being studied, 
and the velocity of the air. Since the density of air in the 
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atmosphere varies with height above sea level, an arbitrary- 
height of 30,000 meters was assumed for all M values. The 
density and temperature of air at this altitude was taken 
from the International Civil Aviation Organization (ICAO) 
standard atmosphere tables [Ref 6]. The calculations may be 
used for any altitude by inserting the appropriate density 
and temperature. 

D. EXERGY CALCULATION 

From equation (1), and neglecting all terms on the right 
except the entropy term, the total exergy decrease through 
the shock wave is proportional to the increase in entropy 
across the shock. The entropy change across a shock is given 
by [Ref. 7:p. 10} 

As 2 T Mfsln J «-W-l) , 

— = ln(---) - yln(-—--—--) 

c v T + 1 ( 7 -DMjsin 2 6 + 2 ( ' 

where c v is the specific heat at constant volume for air, 7 
is the ratio of specific heats, M is the freestream Mach 
number, and 6 is the local shock wave angle with reference 
to the freestream direction. By using the method of Moeckel, 
modified by Love, to determine the approximate shock shape, 
the entropy increase across the total shock was calculated 
by equation 5 and numerical integration using cumulative 
applications of Simpson's Rule. The number of integrations 
varied with the shock strength and shape running from a low 
of 50 for Mach*1.1 to a high of 2513 for Mach=10. 
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The entropy calculations were performed by the Fortran 
77 program contained in Appendix A. For the assumed body 
shape, this program calculates the approximate shock wave 
shape and the total rate of entropy change across the shock 
for a designated control volume and designated Mach numbers 
between 1.1 and 10. The calculations are carried to three 
decimal places with the knowledge that assumptions and 
round-off error may have affected the accuracy. The results 
are summarized in Table 2. 


TABLE 2: ENTROPY CHANGE ACROSS SHOCK 


Mach Number 

Rate of Entropy Change 
Across Shock (KJ/sec} 

l.l 

. . 9 

2.0 

50.803 

3.0 

301.469 

4.0 

836.941 

5.0 

1662.150 

6.0 

2764.771 

7.0 

4131.494 

8.0 

5759.238 

9.0 

7641.229 

10.0 

9774.107 


As can be seen from equation 5, the entropy rise across 
the shock increases as the shock wave angle (6 ), and the 
freestream Mach number (M x ) increase. For a detached shock 
wave, the shock wave angle is greatest immediately ahead of 
the body where it is approximately normal to the freestream 
and decreases as it curves around the body until it 
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eventually asymptotes the freestream. This indicates that 
for a given body shape and Mach number, the greatest exergy 
loss occurs in the vicinity of the vertex of the shock wave. 
This is shown graphically in Figure 3. The vertical axis is 
shown as entropy change over the specific gas constant <R'), 
and the horizontal axis is shown as the angle phi ( <f>) , which 
is 90 - 6. 



Figure 3: Entropy Rise Across Shock Wave. 


E. VERIFICATION OF THE COMPUTER MODEL 

Comparisons between theoretical shock location using the 
program of Appendix A and actual shock location as 
determined by Love for Mach numbers of 1.94 and 6.8 for a 
combination sphere and cylinder body shape are shown in 
Figures 4 and 5 respectively. These comparisons show good 
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agreement between calculated shock shape and shock location. 
The data points corresponding to experimental results were 
taken from charts [Ref. 4:pp 46,51] and therefore, precise 
quantitative evaluations regarding goodness of fit are not 
possible. However, the experimental data appear to be within 
five percent of the predicted values of the shock profile. 



Shock Shape. Mach 1.94 



Shock Shape. Mach 6.8 
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Another reality check for the computer model was made by 
calculating the mass flow of air through the control volume 
defined by the shock wave at l.lsMslO at an altitude of 
30,000 meters. Since the vehicle is supersonic, there is no 
deflection of the flow around the body prior to passing 
through the shock. The validity of the shock shape generated 
by the program is good as shown in Figures 4 and 5 above. 

The shock wave may be thought of as a blunt cone. If the 
cone is terminated at any point, the distance from the 
centerline to the shock is the radius of the base. The mass 
flow through the area of the base can be easily calculated. 
This was done and the result was compared with the mass flow 
calculated by the program using the same Simpson's rule 
integration technique that was used to calculate the entropy 
rise across the shock. The results of these two calculations 
differed by less than three per cent for all Mach numbers. 
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III. EXERGY DECREASE IN THE BOUNDARY LAYER 


If a control volume is taken around a cylindrical body 
moving through a gas, there will be a change in the momentum 
of the gas taken upstream of the body from that taken 
downstream of the body. This difference in momentum is a 
function of the drag on the body, and is given by [Ref. 

8 :pp. 23-24] 


D = / pu(y-u) dA (6) 

J A 

where the integral is taken over the end surface of the 
cylinder at a large distance downstream, U is the velocity 
upstream, and u is the velocity at an arbitrary point in the 
wake. Letting 


U-u - -f dz) 

J o ay dz 

According to Crocco's theorem ^ef. 9] 


(7) 


ay dz 


1 

u 


(pdy 

dy 


ds 

dz 


dz) 


(8) 


Substituting equations 7 and 8 into equation 6, we get 


D= f pudAf — (-f^ dy* dz) (9) 

Ja Jo u dy dz 

Replacing pudA {the mass flow through the cross - sectional 
area dA) by dm and T/u by T m /U, which is justified if second 
order terras are neglected [Ref. 8:p.24], in equation 9 gives 
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D = I=[ dm(s-sj ( 10 ) 

U J A 

Taking the integral over the area of the wake, equation 10 
can be written as 

DU = T m —— ( 11 ) 

sec 

This is an application of the theory of Oswatitsch [Ref. 

10] . The freestream temperature and velocity are known or 
can be calculated using the ICAO tables [Ref. 6] for a given 
altitude and Mach number. If the drag can be calculated, 
then the entropy increase can be calculated using equation 
11. To calculate the profile drag of the body, the flow 
parameters at the boundary layer edge and the boundary layer 
shape are needed. 

The theory that the flow surrounding an immersed body 
could be divided into two parts was first proposed by L. 
Prandtl in 1904. Prandtl stated that the viscosity in the 
major part of the flow field has a negligible effect on its 
motion and consequently in most cases could be ignored. The 
other part of the flow field consisted of a narrow region 
next to the body where the viscous forces are not 
negligible. [Ref.11] 

This narrow region near the body is known as the 
boundary layer. In the boundary layer the velocity of the 
fluid increases rapidly from zero at the surface to a value 
which then changes slowly with increasing distance from the 
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body. At an infinite distance from the body, the boundary 
layer and the freestream merge smoothly together. For 
practical applications, the boundary layer conditions 
approach the freestream condition very rapidly within a 
relatively small distance from the surface as shown in 
Figure 5. 



If the freestream velocity is taken as U, then the edge of 
the boundary layer may be designated as that distance from 
the body where the velocity of the stream is .995U or u e . 

Boundary layers may be either laminar or turbulent. In 
laminar flow, the flow within the boundary layer is 
characterized by smooth streamlines approximately parallel 
to the surface and there is little mixing. In the turbulent 
boundary layer the flow experiences irregular fluctuations 
in velocity magnitude and direction superimposed on the mean 
flow, which remains roughly parallel to the surface 
[Ref. 12:p.11]. 

The change from laminar to turbulent flow is referred to 
as transition. Transition may be caused by a number of 
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factors including disturbances that occur in the freestream, 
surface roughness, changes in the pressure distribution on 
the surface of the body, or an increase in Reynolds number. 
Reynolds number is a ratio of the inertia to viscous forces 
in a fluid and is given by the following equation: 

Re = (12) 

V- 

where U is the velocity, p is the density, 1 is the body 
dimension, and p is the dynamic coefficient of viscosity of 
the fluid. 

The presence of shear forces in the boundary layer due 
to friction and the viscous nature of the fluid result in 
heat being generated. The heat is transferred to the body 
and to the fluid. The viscous dissipation as heat in the 
boundary layer results in an increase in entropy and 
consequently a decrease in exergy. At hypersonic velocities, 
the intense viscous dissipation may result in temperatures 
high enough to cause dissociation and ionization of the 
surrounding air [Ref. 5:p. 228]. This phenomenon and its 
effects will not be discussed in this paper. However, the 
error introduced by ignoring these effects can be 
significant. Interpolating from Anderson [Ref. 5:p. 511], 
the error, at M=10 and at 30,000 meters, in the ratio T 2 /T x 
across a normal shock is approximately 24 per cent. 
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A. BOUNDARY LAYER CALCULATIONS 


An exact determination of the flow conditions in and 
surrounding the boundary layer for a body of arbitrary shape 
is possible by solving the three dimensional Navier-Stokes 
equations. These equations do not have a closed form 
solution, but require elaborate computer methods to be 
solved. However, if we specify that the boundary layer 
thickness (6) is much less than the length of the body, and 
that the Reynolds number is large (>10 5 ), the Navier-Stokes 
equations can be reduced to the boundary layer equations: 


Continuity. . !<££>. - 0 

ox ay 


(13) 


x Momentum: pu-f^+pv-f^ = - ^ ( } (14) 

ox ay dx dy dy 


y Momentum: 


|E = o 

dy 


(15) 


Enetgy. p U |* .pv|* (*|I) + .„<«“) 

dx dy dy dy dx dy 


(16) 


where x is the direction parallel to the surface, y is the 
direction normal to the surface, u is the velocity component 
in the x direction, v is the velocity component in the y 
direction, p is the density, p e is the pressure at the outer 
edge of the boundary layer, p is the viscosity coefficient 
of the fluid, h is the enthalpy, and k is the thermal 
conductivity of the fluid. 
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Like the Navier-Stokes equations, the boundary layer 
equations are nonlinear. They are simpler and easier to 
solve, but still require computer methods to obtain exact 
solutions in most cases. Solution methods for the boundary 
layer equations are presented by Young [Ref. 12] and 
Schlichting [Ref. 14]. 

1. Approximate Methods 

The purpose of this thesis was not to determine an 
exact solution to the boundary layer equations, but to show 
that an exergy analysis can be conducted on the boundary 
layer and to provide the foundation for an exergy based 
design tool and or methodology. This concept can be 
demonstrated by making several approximations that will 
provide rough estimates of the required variables without 
the extensive numerical and computer calculations required 
to obtain a more exact solution of the boundary layer 
equations. The necessary calculations are developed in the 
following paragraphs. 

a. Stagnation Point Calculation 

As stated in the preceding chapter, the flow 
field surrounding a blunt body traveling at supersonic 
velocity will consist of both subsonic and supersonic flow 
divided by the sonic line (see Figure 1). For a body at zero 
angle of attack, the flow precisely along the body axis will 
pass through the shock wave, which at that point is normal 
or very nearly normal to the freestream direction, slow, and 
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3top at the stagnation point. On a blunt body, the boundary 
layer has a finite thickness at the stagnation point. The 
boundary layer thickness at the stagnation point of a two 
dimensional body is given by 


6 ( 0 ) 



0.075U 

f>V 


(17) 


where 6(0) is the boundary layer thickness at the stagnation 
point, and V is the magnitude of the velocity gradient at 
the stagnation point. [Ref. 15:pp. 208-209] 

The velocity gradient is determined by the 
method of Anderson [Ref. 5:pp. 255-256]. It can be shown 
using Euler's equation applied at the boundary layer that, 


du a = l dp a 
dx ~ p a u a dx 


(18) 


where the subscript e indicates conditions at the edge of 
the boundary layer, u is the velocity component in the x 
direction which is parallel to the body, and p is the 
pressure. Assuming a modified newtonian pressure 
distribution over the surface, the pressure coefficient is 
given by [Ref. 5:p. 53] 

C p = C^sin 2 ©, (19) 

where $ x is the angle between the tangent to the surface and 
the freestream direction, and is given in terms of y 

and the freestream Mach number (M x ) as 
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(7 + 1 ) 2 Ml 


( 20 ) 


( [ 


7 m \ 47ftfJ-2(7-l) 


7 + 1 


Defining $ as the angle between the normal to the surface 
and the freestream, then equation 19 can be written as 


c p = c *w C0s2 $ 


( 21 ) 


Defining C p as 


r _ P*-P i 

P rr 


( 22 ) 


where the subscript 1 denotes the freestream condition, then 
equation 21 can be written as: 


Pa — = C cos 2 <|> 

q l 


(23) 


where q x is the dynamic pressure of the freestream. Equation 
23 can be written as 


Pe 9 


(24) 


Differentiating equation 24 with respect to x, we obtain 




Substituting equation 25 into equation 18, we have 
du c 


dx 


— ft** — cos<|>sin<t> — 


P e “ e 


dx 


(25) 


( 26 ) 


Consider the stagnation region, as shown in Figure 7. 
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Figure 7: Stagnation Region of a Blunt Body [Ref.5:p. 
251] . 

In this region, let Ax be a small increment of surface 
distance above the stagnation point, corresponding to small 
changes in <f> of A$. Since the flow velocity is so low in the 
stagnation region, we can assume [Ref. 5:p. 252] 

u* = (-$ 2 > Ax (27) 

* dx s 

In the stagnation region the following approximations can be 
made if <f> is small: 


COS<|> = 1 

(28a) 

sin<t> =<J)sA4> = -^ 

fv 

(28b) 

n 

Sl-S 

(28c) 


where R is the radius of the nose of the body at the 
stagnation point. 


24 









(31) 

This value for the velocity gradient is then substituted 
into the value for V in equation 17. 

The velocity, temperature, density, Mach number, 
and pressure immediately behind the shock wave can be 
determined from the normal and oblique shock relations 
contained in Anderson [Ref. 2:pp.67-68,106]. The density 
behind the shock along the stagnation streamline is used in 
equation 17 for the value of p. 

Equation 17 gives the boundary layer thickness 
at the stagnation point of a two dimensional body. The 
boundary layer thickness for an axisymmetric body will be 
thinner than for the two dimensional case [Ref. 5:p. 254]. 
While the difference between the axisymmetric and two 
dimensional cases is noted, equation 17 does provide an 
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approximation of the boundary layer thickness for the 
axisymmetric case at the stagnation point. 

b. Boundary Layer Edge Conditiona 

The pressure distribution at the surface of the 
body is approximated by the modified newtonian method [Ref. 
5:pp. 53-56]. Although the results provided by this method 
are not very accurate at low supersonic Mach numbers, the 
results improve quickly as Mach number increases [Ref.5:pp. 
45-56]. The pressure coefficient at the surface of the body 
using the modified newtonian method is given by equation 19. 
Recalling from equation 22 the definition for C p/ the 
pressure at the boundary layer edge over the nose of the 
body is approximated by 

p e = C p<gu q 1 co3 2 ff *p 1 (32) 

Figure 8 shows the surface pressure distribution over the 
nose of the hypothetical vehicle used in this analysis for 
selected Mach numbers. The nose radius (R) in Figure 8 is 
1.5 ft (.4572 m) and the pressure distribution is shown from 
the stagnation point through an angle of 68.2 degrees, where 
the rounded nose intersects the cone of the body. 

The total temperature behind the shock wave 
remains constant everywhere behind the shock. Also, we 
assume that the total pressure is constant along the 
streamline that enters through the nearly normal portion of 
the shock and flows along the boundary layer edge. These 
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values can be calculated and used as a reference to obtain 
local values of temperature and pressure. 



Figure 8: Pressure Distribution Over a Blunt Nose 
Cylindrical Cone Body. 


The total pressure in the freestream is given by 

P 0i = + (33) 

where P 01 is the stagnation pressure in the freestream, P is 
the static pressure of the freestream. For these analyzes 
all freestream atmospheric conditions were taken from 
Reference [6] at an altitude of 30,000 meters. The ratio of 
stagnation pressures across the shock is 
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*0. 

p 0, 


e 


s , - s, 

R' 


(34) 


where R' is the specific gas constant, and s 2 -s 1 is given by 


s 2 -s 1 = CplnKl + il -1)) 


-Rln[l + JiL (Ml -1) ] 

“V • L 


2 ♦ (y-l)M{ 
(y + i )M?> 


) 3 


(35) 


where c p is the specific heat at constant pressure of the 
gas. The stagnation pressure behind the shock is then 
obtained by taking the product of equations 33 and 34 


Pn = 



(36) 


The stagnation temperature in the freestream is 
calculated by 

T 0i = T(l + IZl/f?) (37) 

where T is the freestream temperature. The total temperature 
across a normal shock is constant, therefore 


T 0 , = t o 2 ( 38 ) 

If we substitute T 02 for T 01 in equation 37, and with a 
known Mach number, the local temperature can be determined. 

The sonic point on the body nose is located with 
the same method used in determining the shock shape, and is 
described in Appendix A. Since the Mach number at the sonic 
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line is one, Che local pressure at the sonic point can be 
obtained by substituting the value of P 02 , from equation 36, 
for P 01 in equation 33. Once the sonic point is located, the 
angle between the body centerline and the sonic point can be 
easily calculated. This angle will be designated the sonic 
point angle (SPA) for calculation purposes. The total angle 
that the flow makes as it rounds the nose is 68.2 degrees. 
Subtracting the SPA from 68.2 gives the remaining angle, 
designated A 1# that the flow must turn through in traveling 
from the sonic line to the cone portion of the body. 

Treating angle A x as a Prandtl-Meyer expansion angle, the 
pressure, Mach number, and temperature of the flow 
downstream of the expansion can be calculated. 

Before proceeding further, it is necessary to 
introduce the Prandtl-Meyer function [Ref. 2:p. 134], which 
is given the symbol v 


v (M) 


111 tan" 1 lZl(M 2 - 1) -tan"VM 2 -l 
7-1 y +1 


(39) 


As can be seen from equation 39, the Prandtl-Meyer function 
is a function of the Mach number and y. Let v (M-^) be the 
Prandtl-Meyer function before the expansion, and MM 2 ) be 
the Prandtl-Meyer function after the expansion. Then the 
angle A 1 that the flow turns through is 
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At the sonic line the Mach number is one, and 
therefore i»(M 1 > in equation 40 is zero and A r MM 2 ) . Putting 
the known value of A x into equation 39 for v (M), the Mach 
number after the expansion can be determined. This process 
is repeated at the point wb"?re the cone joins the body. In 
this fashion the flow Mach number along the length of the 
cylindrical portion of the vehicle can be found, with the 
Mach number known and the previously calculated values of 
total pressure and total temperature, the local pressure and 
temperature along the cylindrical portion of the body can be 
determined using equations 33 and 37. These values will be 
taken as the conditions at the edge of the boundary layer, 
c. Boundary Layer ThlckneaB 

The boundary layer is formed as a result of the 
no slip condition that exists at the wall of a body immersed 
in a flow (see Figure 6). As the flow progresses along the 
length of the body, the thickness (5) of the layer of 
retarded flow increases as more and more fluid becomes 
affected [Ref. 14:p. 25]. The thickness of the boundary 
layer on a flat plate at zero incidence can be calculated 
using the method developed by Blasius [Ref. 16]. For laminar 
flow the boundary layer thickness was found to be 



where x is the distance along the plate, and Re x is the 
Reynolds number (see equation 12) as a function of distance 
along the plate. For a turbulent boundary layer the 
thickness is given by [Ref. 12:p. 13] 


0.37x 
(Re x 5 1/5 


(42) 


As previously stated, the boundary layer on an 
axisymmetric body is thinner than on a flat plate. However, 
as a first approximation, the boundary layer over a flat 
plate with a length equal to the distance traveled by the 
flow from the stagnation point to the tail end of the body 
will be used to approximate the boundary layer thickness. 
This approximation is reasonable when 6<<R as it is in this 
case. 

Since there is a dramatic difference in the 
growth of laminar vs turbulent boundary layers, it is 
necessary to establish the point where transition occurs. In 
actuality transition may occur over a finite distance making 
prediction of a point of transition extremely difficult. 
Research to determine exactly how transition occurs and 
methods to predict its occurrence are ongoing. A formula 
developed by Anderson [Ref. 5:p 280] to predict the point of 
transition will be used in this analysis. 
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( 43 ) 


log 10 (J2e T ) = 6.421 exp [1.209x10 ~ 4 Mg • 641 ) 

where Re T is the Reynolds number at which transition occurs, 
and M e is the Mach number at the boundary layer edge. It 
should be recognized that the data used to develop this 
equation was unknown to the author and that equation 42 was 
being used as an approximation. 

B. BASE PRESSURE DETERMINATION 

If the assumption is made that the launch vehicle is a 
projectile with no exhaust plume, the pressure behind the 
body in the wake is assumed to be the freestream static 
pressure. This assumption is supported by Anderson [Ref. 

5:pp 117-138] who shows that for a cylindrical body with 
fineness ratio greater than six, the pressure asymptotically 
approaches freestream pressure. The assumption of no exhaust 
plume allows calculation of the entropy change as a result 
of the boundary layer on the body, but ignores the effects 
of boundary layer/plume interaction. 

With the pressure in the wake being approximately equal 
to static pressure and the stagnation pressure behind the 
shock, from equation 36, the Mach number in the wake can be 
determined using equation 33. The difference between the 
Mach number in the wake and in the flowfield can be used to 
calculate the difference between the velocities (U-u) needed 
to solve equation 6. These values are summarized in Table 3. 
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C. WAKE PROFILE 


The wake boundary is established from the body plus the 
boundary layer thickness at the base of the body. Traveling 
downstream for several body diameters and assuming the same 
rate of growth of the boundary layer, the distance from the 
centerline to the point where the wake velocity approaches 
the freestream velocity can be determined. This is shown in 
Figure 9. For calculation purposes, a distance of five 
diameters downstream of the base was arbitrarily chosen. 

TABLE 3: COMPARISON OF FREESTREAM AND WAKE MACH NUMBER 


Freestream Mach (M^) 



1.1 

1.0991 

. , 3 92 

2.0 

1.7879 

.8939 

3.0 

2.2769 

.7590 

4.0 

2.6351 

.6588 

5.0 

2.9218 

.5844 

6.0 

3.1630 

.5272 

7.0 

3.3725 

.4818 

o 

CD 

3.5586 

.4448 

9.0 

3.7265 

.4141 

10.0 

3.8800 

.3880 



Figure 9: Illustration of Wake Profile. 
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The variation of u in the wake from the centerline to 
the point where u»U is given by 

u = u w * x(i--^r) y (44) 

where U w is the minimum wake velocity, U x is the freestream 
velocity, and x is a scaling variable where Osxsl. The 
distance from the centerline to the point where u~U is given 
by 

y = (R b + 8 C ) n (45) 

where R b is the radius of the base of the body and 6 t is the 
thickness of the boundary layer extended to the point in the 
wake where the profile is being calculated. 

D. DRAG DETERMINATION 

The drag of the body is calculated using equation 6 
which is repeated here: 

D = f pu(U-u)dA (46) 

J A 

Since the velocity in the wake varies along the centerline, 
designated the x axis, and with the distance from the 
centerline, designated the y axis, equation 46 becomes the 
following double integral: 

D - f y f 1 2nrpu(U-u)dr (47) 

where u and y are given by equations 44, and 45 
respectively. The density in equation 47 is the density at 
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the location in the wake where the profile is determined. 

For these calculations this value was approximated using the 
perfect gas equation. The pressure is taken as freestream 
and the temperature is that calculated for the boundary- 
layer temperature just forward of the base of the body 
ignoring plume effects. 

E. ENTROPY CALCULATION 

With the drag determined, the entropy flux increase as a 
result of the boundary layer can now be calculated using 
equation 11 which is repeated here 

DU = 71 —— ( 48 ) 

1 sec 

where U and T x are the freestream velocity and temperature 
respectively, and D is the drag. 

For the body shape used in this analysis and with the 
approximations used to calculate the required flowfield 
parameters, the entropy increase is presented in Table 4. 

As can be seen from Table 4, the entropy increase in the 
boundary layer starts to decrease above M=7. This occurs due 
to the boundary layer becoming thinner at higher Ma li 
numbers. The thinner boundary layer is a result of the 
transition point moving aft on the body. Initially the 
transition point moves forward toward the nose of the body. 
Between 6sMs7, the transition point stops moving toward the 
nose and starts moving back toward the tail. This is 
occurring as a result of the laminar to turbulent 
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TABLE 4: ENTROPY INCREASE DUE TO THE BOUNDARY LAYER 


Mach number 

Rate of Entropy increase 
(KJ/sec) 

1.1 

2918 

2 

21230 

3 

48316 

4 

70635 

5 

84925 

6 

92074 

7 

92671 

8 

87442 

9 

77123 

10 

62372 


transition point as calculated by equation 43. The 
transition Reynolds number increases exponentially with 
increasing Mach number. For a given transition Reynolds 
number, the point of transition is obtained by substituting 
Re T for Re and x for 1 in equation 12, and rearranging 


x 


Re ipPL 
~PV~ 


(49) 


As Re T increases exponentially, the numerator in equation 49 
grows faster than the denominator, as shown in Figure 10. 

The result of this is that the transition point moves toward 
the tail of the body and the boundary layer remains laminar 
longer. Since entropy is produced at a higher rate in 
turbulent boundary layers, the net result is that less 
entropy production occurs. 
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Mach Number. 


This phenomenon is a result of using equation 43 to 
predict the transition point and may not be an accurate 
representation of actual flow. As previously stated, the 
data that equation 43 was developed from was unknown to the 
author, and the conditions of this study may be outside the 
range for which equation 43 is valid. 


37 
















IV. CONCLUSIONS 


Recall from the exergy equation (equation 1) that an 
increase in entropy results in a corresponding decrease in 
exergy. A primary advantage of exergy analysis is the 
identification of those processes or areas that have the 
highest exergy loss and hence the greatest potential for 
improved efficiency. A method to quantify the exergy 
decrease resulting from entropy production across the shock 
wave and in the boundary layer, ignoring high temperature 
effects, has been presented. 

It has been demonstrated that the exergy decrease due to 
the production of entropy across the shock wave and in the 
boundary layer can be calculated for a typical space launch 
vehicle. This information can be used to highlight those 
areas which have the greatest potential for improved 
efficiency. 

For the conditions of this study and within the accuracy 
of the approximations, it is clearly evident that entropy 
production :.s much greater in the boundary layer than across 
the shock wave for a given Mach number. This indicates that, 
all other factors being equal, design efforts should 
concentrate on reducing the entropy production in the 
boundary layer. 
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V. RECOMMENDATIONS POR FURTHER RESEARCH 


This thesis applies the method of exergy analysis to two 
areas of interest for space launch vehicles. In doing so, it 
provides an initial step in the development of an exergy 
methodology tool for design purposes. 

Following the work of Czysz [Ref. 21], application of 
this methodology to launch vehicle propulsion systems would 
be beneficial. 

Alternately, variations in launch vehicle design could 
be examined with the goal of finding the optimum design to 
minimize entropy production across the shock wave and within 
the boundary layer. 

It would also be useful to determine the magnitude of 
the error for the approximations used in this study. This 
could be accomplished by comparing the results with those 
obtained from exact flowfield calculations using CFD 
methods. Also the effects of the boundary layer/plume 
interactions on entropy production could be investigated. 
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APPENDIX A 
THE SHOCK WAVE MODEL 


This program computes the entropy increase across the 
shock wave of a cylindrical space launch vehicle with a 
blunt cone nose. The program is written for Mach numbers of 
1.1 to 10, and for atmospheric conditions at an altitude of 
30,000 meters. This program is written in Fortran 77 code. 
Comment statements are preceded by a "C". 


C Variable declarations 

M = Freestream Mach number 
BETA » Cotangent of the Mach number 

THETA = The shock wave angle for sonic flow behind the 
shock 

DELTA - The maximum stream deflection angle that can 
occur without shock separation 
NU = Angle between the sonic line and a line normal to 
the freestream direction 
GAMMA - Ration of specific heats (1.4 for air) 

SIGMA * Angle of shock wave at the center of mass 
streamline 

A = Area ratio for isentropic flow 
P = Stagnation pressure ratio 

RATIO » Ratio of the distance from the x-axis to the 

shock sonic point to the distance from the x-axis 
to the body sonic point 

YSB * Distance from the centerline of the body to the 
sonic point on the body 

X0 - Distance from foremost point of detached shock to 
intercept of its asymptote on the x-axis 
Y » Coordinate perpendicular to freestream flow 
X a Coordinate in the direction of the freestream 
I,J,Z ■ Increment counters 
PI a Radian equivalent of 180 degrees 

DIFF - The difference between the shock wave angle and 
the asymptote angle of the hyperbola that is 
defined by the shock wave 

CANGLE » The constant angle between the asymptote of 
the hyperbola and the freestream direction 
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C SANGLE = The angle between the freestream direction and 
C the tangent to the shock wave 

C ANGLED = The same angle as "SANGLE" but in degrees 
C instead of radians 

C TOTALS = The total change in entropy across the shock 
C wave 

START PROGRAM 
INTEGER I, J, Z 
PARAMETER (Z-10) 

REAL M(Z) , BETA(Z), THETA(Z) , DELTA(Z), 

+ NU(Z), PI, YSB(Z), SIGMA(Z) , A{Z), P(Z), 

+ RATIO(Z), X0(Z), X(7500), Y{7500), DIFF, B(Z), 

+ CANGLE{Z), SANGLE(7500), ANGLED(7500), TOTALS 

PI=3.14159 
M (1)=1.1 
M {2 ) =2 

BETA (1) =SQRT (M (1) **2.0-1.0) 

BETA(2)=SQRT(M(2)**2.0-1.0) 

DO 05 1=3,Z 

M(I)=M(I-1)+1.0 

BETA(I)=SQRT(M(I)**2.0-1.0) 

05 CONTINUE 

CALL SHOCK(Z,M,THETA) 

CALL MAX(Z,M,DELTA) 

CALL AVG(Z,M,NU) 

CALL DISYSB(Z,DELTA,YSB) 

CALL CENTER(Z,THETA,BETA,SIGMA) 

CALL AREA(Z,M,A) 

CALL PRESS(Z,SIGMA,M,P) 

CALL SPRATIO(Z,P,A,NU,RATIO) 

CALL DISXO(Z,YSB,BETA,RATIO,THETA,X0) 

C The following do loop, for each Mach number, calculates 
C the angle of the asymptote of the hyperbola that 

C defines the shock. Then the point on the x and y 

C axis where the difference between the shock wave 
C angle and the asymptote angle is less than 1 per 

C cent is found. This defines the limit of the shock wave 

C for purposes of calculating the entropy change across 
C the shock. The change in entropy is then calculated. 

DO 15 J=l,Z 

B(J)=X0{J)/BETA(J) 

CANGLE (J) = ATAN {B (J) /X0 (J) ) 

CALL SHAPE(X0(J),BETA(J),X,Y,) 

CALL SWANGLE(X0(J),BETA(J),Y,SANGLE,ANGLED) 

DO 10 1=2,7500 

IF (SANGLE(I).EQ.0) GO TO 10 
DIFF=(SANGLE(I)-CANGLE(J))/CANGLE(J) 

IF(DIFF.GT.0.01) GO TO 10 
IF(DIFF.LE.0.01) THEN 

CALL TOTAL(M(J),X0(J),BETA(J),X(I),TOTALS) 

GO TO 15 
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END IF 


10 CONTINUE 
15 CONTINUE 
STOP 
END 


C Calculate the shock wave angle for sonic flow behind 
C the shock 

INTEGER I,Z 

REAL M(Z), THETA(Z), GAMMA 
GAMMA=1.4 
DO 05 1=1,Z 

THETA(I)=ASIN(SQRT((1.0/4.Q*GAMMA*M(I)**2.0)) 

+ *((GAMMA+1.0)*M(I)**2.0-(3.0- GAMMA) 

+ +(SQRT((GAMMA+1.0)*((GAMMA+1.0)*M(I)**4.0-2.0 

+ *(3.0-GAMMA)*M(I)**2.0+(GAMMA+9.0))))))) 

05 CONTINUE 
RETURN 
END 

C Calculate the maximum deflection angle, in radians, 

C for a given freestream Mach number in air. For this 

C calculation, the maximum deflection angle without shock 

C detachment was calculated from NACA report 1135, chart 
C 5. An equation for the line of maximum deflection was 

C obtained by performing a curve fit of the data 

C presented in the chart. 

SUBROUTINE MAX(Z,M,DELTA) 

C NLOG = Natural log of Mach number between l.l and 10 

C TVAR = Natural log of shock wave detachment angle in 

C degrees 

INTEGER I,Z 

REAL M(Z), NLOG(10), TVAR(10), DELTA{Z) 

DO 05 1=1,Z 

NLOG(I)= LOG(M(I)) 

TVAR(I)=3.80568+0.41631*L0G(NLOG(I)) 

DELTA(I)=(EXP(TVAR(I)))/57.29578 
05 CONTINUE 
RETURN 
END 


C Calculate the angle between the sonic line and a line 
C normal to the freestream direction 
SUBROUTINE AVG(Z,M,NU) 

INTEGER I,Z 
REAL M(Z), NU(Z) 

DO 05 1=1,Z 

NU(I)=(3.9901+56.9298*LOG(M(I)))/57.29578 
05 CONTINUE 
RETURN 
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END 

C Calculate the distance from the body centerline to the 
C sonic point on the body 

SUBROUTINE DISYSB(Z,DELTA,YSB) 

C DISTl=Distance along a line drawn tangent to the body 
C sonic point from the x-axis to the sonic point on 

C the body 

C R=Radius of the body apex (in feet) 

C ANG3=The angle between the x-axis and a line 
C perpendicular to DIST1 at the body sonic point and 

C extending from the sonic point to the x-axis 

INTEGER I,Z 

REAL PI, YSB(Z), DISTl(lO), R, ANG3(10), DELTA(Z) 

PI=3.14159 
R-1.5 

DO 05 1=1,Z 

ANG3(I) = PI-(PI/2.0)-DELTA(I) 

DIST1(I)=R*((SIN(ANG(I))/SIN(DELTA(I)))) 

YSB(I)=DIST1(I)* SIN(DELTA(I)) 

05 CONTINUE 
RETURN 
END 

C Calculate the angle of the shock wave at the center 
C of mass streamline between the shock wave vertex 

C and the sonic point of the shock wave 

SUBROUTINE CENTER(Z,THETA,BETA,SIGMA) 

INTEGER I, Z 

REAL THETA(Z), BETA(Z), SIGMA(Z) 

DO 05 1=1,Z 

SIGMA(I)=ATAN(((TAN(THETA(I)))**2.0*(3.0/2.0)**2.0 
+ -(1.0/BETA(I))**2.0)*((3.0/2.0)**2.0-1.0))**0.5) 

05 CONTINUE 
RETURN 
END 

C Calculate the area ratio for isentropic flow 
SUBROUTINE AREA(Z,M,A) 

INTEGER I,Z 

REAL GAMMA, M(Z), A(Z) 

GAMMA=1.4 
DO 05 1=1,Z 

A(I)=((GAMMA+1.0)/2.0)**((GAMMA+1.0)/(2.0 
+ *GAMMA-1.0)))*M(I)*(1.0+((GAMMA-1.0)/2.0) 

+ *M(I)**2.0)**(-((GAMMA+1.0)/(2.0*(GAMMA-1.0)))) 

05 CONTINUE 
RETURN 
END 
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C Calculate the ratio between stagnation pressure before 

C and after the shock wave for a given shock wave angle 

SUBROUTINE PRESS(Z,SIGMA,M ?) 

INTEGER I, Z 

REAL GAMMA, SIGMA(Z), M(Z), P(Z) 

GAMMA = 1.4 
DO 05 1=1,Z 

PI(I)=1.0/{ ( ((GAMMA+1.Q)*M(I)**2.0* 

+ SIN(SIGMA(I))**2.0)/((GAMMA-1.0)* 

+ M(I)**2.Q*SIN(SIGMA(I))**2.0+2.0))**(GAMMA/ 

+ (GAMMA-1.0))*((GAMMA+1.0)/(2.0*GAMMA*M(I)**2.0* 

+ SIN(SIGMA(I))**2.0-(GAMMA-1.0)))** 

+ (1.0/(GAMMA-1.0))) 

05 CONTINUE 
RETURN 
END 

C Calculate the ratio of the distance from the x-axis to 
C the shock sonic point to the distance from the x- 
C axis to the body sonic point 

SUBROUTINE SPRATIO(Z,P,A,NU,RATIO) 

INTEGER I,Z 

REAL P(Z), A(Zj , NU(Z), RATIO 
DO 05 1=1,Z 

RATIO(I)= 1.0/((1.0-P(I)*A(I)*COS(NU(I)))**0.5> 

05 CONTINUE 
RETURN 
END 

C Calculate the distance from the foremost point of the 

C detached shock wave to the intercept of its asymptote 

C on the x-axis 

SUBROUTINE DISX0(Z,YSB,BETA,RATIO,THETA,X0) 

INTEGER I, Z 

REAL YSB(Z), BETA(Z), RATIO(Z), THETA{Z), X0(Z) 

DO 05 1=1,Z 

X0(I)=YSB(I)*BETA(I)*RATIO(I)*((BETA(I))**2.0* 

+ (TAN(THETA(I)))**2.0-1.0)**0.5 

05 CONTINUE 
RETURN 
END 

C Calculate the x and y coordinates for any point on the 
C shock wave 

SUBROUTINE SHAPE(X0,BETA,X,Y) 

INTEGER I 

REAL X(7500), Y(7500), BETA, X0 
X(l)=0.1 
Y(1)=0 

DO 05 1=2,7500 

X(I)«X(I-1)+0.2 
IF(X(I).LT.X0) THEN 




Y(I)=0 

ELSEIF (X(I).GE.XO) THEN 
Y{I)=SQRT(X(I)**2.0-XO **2.0)/BETA 
END IF 
05 CONTINUE 
RETURN 
END 

C Calculate the angle between the freestream and a line 

C drawn tangent to the shock wave at any point 

SUBROUTINE SWANGLE(XO,BETA,Y,SANGLE,ANGLED) 

INTEGER I 

REAL XO, BETA, Y(7500), SANGLE(7500), PI, ANGLED(7500) 
PI-3.14159 
DO 05 1=1,7500 

IF (Y(I).EQ.O) THEN 
SANGLE(I)=0 

ELSEIF (Y(I).GT.O) THEN 

SANGLE(I)=ATAN((SQRT(X0**2.0+BETA**2.C*Y(I)**2.0) ) 
+ /(BETA**2.0*Y{I) ) ) 

ANGLED(I)=SANGLE(I)*180/PI 
END IF 
05 CONTINUE 
RETURN 
END 

Calculate the total change in entropy across the shock 
wave. All calculations are in SI units 
RO=The density of air. For these calculations this 

value is taken from the ICAO standard atmosphere 
tables at an altitude of 30,000 meters. 

TEMP=The temperature in degrees Kelvin of the 

atmosphere at 30,000 meters altitude, taken from 
the ICAO standard atmosphere tables 
GAMMA=Ratio of specific heats 

SGC=Specific gas constant. Assuming an ideal gas, for 
air this value is 287 Joules/Kg-K 
SSA=Local speed of sound in air 

VEL=Velocity, in meters per second, of the freestream 
SHCV=Specific heat of air at constant volume 
PI=Radian equivalent of 180 degrees 

LLINT=Lower limit of integration. The starting point, 
located at the shock wave vertex, for numerical 
integration to calculate the entropy change across 
the shock wave 

ULINT=Upper limit of integration. The ending point for 
numerical integration across the shock wave 
H=The integration step size for application cf 
Simpson's rule 

N=The number of Simpson's rule integrations, and the 
ending point for summation of individual 
integrations 
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CENTERX=The center point, along the x-axis, for each 
Simpson's rule integration 
TSTEP-Specifies the beginning (TSTEP(l)}, center 
(TSTEP ( 2 )), and end (TSTEP{3)) points for 
Simpson's rule application for each individual 
summation 

RADIUS=The distance from the centerline x-axis to the 
shock wave at each TSTEP location 
TANGLE=The shock wave angle at each point of 

intersection of the radius and the shock wave 
DELTAE=Change in entropy per mass flow across the shock 
at each point of intersection of RADIUS with the 
shock wave 

DELTAS=Distance along the shock at each interval of 
integration 

HOLD=Temporary storage variable 
INTEGER I, K 

REAL XO, X, BETA, DELTAE(3), TOTALS, LLINT, ULINT, H, 

+ CENTERX, PI, RADIUS(3), SHCV, GAMMA, TANGLE(3), 

+ TSTEP(3), RO, VEL, SSA, TEMP, M, DELTAS 

R0-.01841 
TEMP=226.509 
GAMMA=1.4 
SGC=287 

SSA=SQRT (GAMMA* SGC* TEMP) 

VEL=M*SSA 

SHCV=717.5 

PI-3.14159 

LLINT=X0*.3048 

ULINT=X*.3048 

HOLD=0 

H=0.5 

N=INT((ULINT-LLINT)/(2.0*H)) 

CENTERX=LLINT+.05 
DO 05 I-LLINT, N 
DO 10 K=1,3 
TSTEP(1)=CENTERX-.05 
TSTEP(2)-CENTERX 
TSTEP(3)=CENTERX+.05 
IF {TSTEP(K).LE.LLINT) THEN 
RADIUS(K)=0 
TANGLE(K)=0 
DELTAE(K)=0 
GO TO 10 
END IF 

RADIUS(K) »(SQRT(TSTEP(K)**2.0-LLINT**2.0))/BETA 
TANGLE(K)=ATAN((SQRT(LLINT**2.0+BETA**2.0*RADIUS(K) 

+ **2.0))/(BETA**2.0*RADIUS(K))) 

DELTAE(K)=SHCV*(LOG((2.0*GAMMA*M**2.0*SIN(TANGLE(K)) 

+ **2.0-(GAMMA-1.0) )./(GAMMA+1.0) )-GAMMA* 

+ LOG(((GAMMA+1.0)*M**2.0*SIN(TANGLE(K))**2.0)/ 

+ ((GAMMA-1.0)*M**2.0*SIN(TANGLE(K))**2.0+2.0))) 
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10 CONTINUE 

DELTAS == (0.1*. 3048) /COS (TANGLE (2) ) 

HOLD=HOLD+(((H/3.0)*(DELTAE(1)+4.0*DELTAE(2) 

+ +DELTAE(3)))*RO*2.0*PI*RADIUS(2)*DELTAS*VEL* 

+ SIN(TANGLE(2))) 

CENTERX=CENTERX+0.1 
05 CONTINUE 

TOTALS=HOLD*TEMP 

RETURN 

END 
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